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ABSTRACT 

We present a technique to explore the radio sky into the nanoJansky regime by employing image stacking 
using the FIRST survey. We begin with a discussion of the non-intuitive relationship between the mean and 
median values of a non-Gaussian distribution in which measurements of the members of the distribution are 
dominated by noise. Following a detailed examination of the systematic effects present in the 20 cm VLA 
snapshot images that comprise FIRST, we demonstrate that image stacking allows us to recover the average 
properties of source populations with flux densities a factor of 30 or more below the rms noise level. With the 
calibration described herein, mean estimates of radio flux density, luminosity, radio loudness, etc. are derivable 
for any undetected source class having arcsecond positional accuracy. 

We demonstrate the utility of this technique by exploring the radio properties of quasars found in the Sloan 
Digital Sky Survey. We compute the mean luminosities and radio-loudness parameters for 41,295 quasars 
in the SDSS DR3 catalog. There is a tight correlation between optical and radio luminosity, with the radio 
luminosity increasing as the 0.85 power of optical luminosity. This implies declining radio-loudness with 
optical luminosity: the most luminous objects (Muv = -28.5) have on average three times lower radio-to-optical 
ratios than the least luminous objects (Muv = -20). There is also a striking correlation between optical color 
and radio loudness: quasars that are either redder or bluer than the norm are brighter radio sources. Quasars 
having g — r ~ 0.8 magnitudes redder than the SDSS composite spectrum are found to have radio-loudness 
ratios that are higher by a factor of 10. We explore the longstanding question of whether a radio-loud/radio- 
quiet dichotomy exists in quasars, finding that optical selection effects probably dominate the distribution 
function of radio loudness, which has at most a modest (~ 20%) inflection between the radio-loud and radio- 
quiet ends of the distribution. We also examine the radio properties of the subsample of quasars with broad 
absorption lines, finding, surprisingly, that BAL quasars have higher mean radio flux densities at all redshifts, 
with the greatest disparity arising in the rare low-ionization BAL subclass. We conclude with examples of other 
problems for which the stacking analysis developed here is likely to be of use. 

Subject headings: surveys — catalogs — radio continuum: general — quasars: general — quasars: absorption 
lines 



1. INTRODUCTION 

"Blank sky" is rarely truly blank. All astronomical imaging 
observations have a sensitivity threshold below which "ob- 
jects" are not detectable. Assuming a reasonably linear de- 
tector response, however, it is not necessarily the case that 
zero photons from discrete sources have been detected at a 
given "blank" spot in an image. If one has reason to believe 
from observations in another wavelength regime that discrete 
emitters are present at a set of well-specified locations, it is 
possible to usefully constrain, or even to detect, the mean flux 
of that set of emitters by stacking their "blank sky" locations. 
The prerequisites for successfully stacking images in this way 
are good astrometry for both the target objects and the survey 
images, and sufficient sky coverage to include a large sample 
of the target class. 

In an early application of this technique, Caillault & 
Helfand (1985) detected the mean X-ray flux from undetected 
G-stars in the Pleiades, using it to constrain the decay of stel- 
lar X-ray emission with age. As higher-resolution X-ray mir- 
rors and detectors have become available over the past two 
decades, X-ray stacking has become a standard analysis tech- 
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nique. Applications have ranged from determining the mean 
X-ray luminosity of object classes in deep X-ray images — 
e.g., normal galaxies (Brandt et al. 2001a), Lyman Break 
galaxies (Brandt et al. 2001b), and radio sources (Georgakakis 
et al. 2003) — to determining the X-ray cluster emission from 
distant clusters in the Rosat All-Sky Survey (Bartelmann & 
White 2003). 

As linear digital detectors have come to dominate optical 
and infrared sky surveys, the stacking technique has been 
widely adopted: e.g., Zibetti et al. (2005) detected intraclus- 
ter light by stacking 683 Sloan Digital Sky Survey (SDSS 6 ) 
clusters, Lin et al. (2004) stacked 2MASS data on cluster 
galaxies, Hogg et al. (1997) stacked Keck IR data to get faint 
galaxy colors, and Minchin et al. (2003) went so far as to stack 
digitized films from the UK Schmidt telescope to comple- 
ment a deep H I survey with the Parkes multi-beam receiver. 
Scaramella et al. (1993) have even stacked cosmological sim- 
ulations in investigating the Sunyaev-Zeldovich effect on the 
cosmic microwave background. 

The radio sky is relatively sparsely populated with sources. 
The deepest large-area sky survey, FIRST, has a surface den- 
sity of only ~ 90 sources deg" 2 at its 20 cm flux density 

6 Funding for the SDSS and SDSS-11 has been provided by the Alfred P. 
Sloan Foundation, the Participating Institutions, the National Science Foun- 
dation, the U.S. Department of Energy, the National Aeronautics and Space 
Administration, the Japanese Monbukagakusho, the Max Planck Society, and 
the Higher Education Funding Council for England. The SDSS Web Site is 
|http : / /www, sdss ■ orq| 
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threshold of 1 .0 mJy. Fluctuation analysis of the deepest radio 
images ever made suggest a source surface density of ~ 15 
arcmin" 2 at ~ 1 /iJy (Windhorst et al. 1993); given that the 
mean angular size of such sources is ~ 2.4", even at these flux 
density levels only ~ 3% of the sky is covered by radio emis- 
sion. Nonetheless, applications of stacking in the radio band 
have been limited. Recently, Serjeant at el. (2004) stacked 
SCUBA data to find the mean submillimeter flux of Spitzer 
24 /xm-selected galaxies. It is with large-area surveys and 
large counterpart catalogs, however, that the stacking tech- 
nique allows one to reach extremely faint flux density levels 
unachievable by direct observations. 

The FIRST survey (Becker, White, & Helfand 1995) is ide- 
ally suited for stacking studies. It contains 811,000 sources 
brighter than ~ 1 mJy over 9030 deg 2 of the northern sky and 
has an angular resolution of 5". Thus, over 99.9% of its five 
billion beam areas represent blank sky. Having written several 
dozen papers on sources detected in the survey, we turn here 
to analyzing the remaining 99.9% of the data. Glikman et al. 
(2004b) presented our initial results of applying radio image 
stacking to the FIRST survey. Wals et al. (2005) applied this 
technique to the 2dF quasar catalog (Croom et al. 2004) using 
FIRST images, producing an estimate of the mean flux den- 
sity for undetected quasars in the range 20-40 /iJy. In this 
paper, we describe the set of detailed tests we have carried out 
to calibrate biases in the FIRST survey's VLA images in order 
to put stacking results on a quantitative basis, and we illustrate 
the technique with several examples. Subsequent papers will 
apply these results to various problems of interest. 

We begin (§2) with a discussion of the median stacking 
procedure that we have adopted, exploring in some detail the 
meaning of, and distinctions between, mean and median val- 
ues in data dominated by noise. We go on to provide a thor- 
ough analysis of the noise characteristics of the FIRST im- 
ages, both by stacking known sub-threshold sources and by 
the use of artificial sources inserted into the survey data (§3). 
We find a non-linear correction to the flux densities derived 
from a stacking analysis, which most likely arises from the 
application of the highly nonlinear 'CLEAN' algorithm to 
these undersampled uv (snapshot) data. We then apply our 
calibrated stacking procedure to the SDSS DR3 quasar sur- 
vey from Schneider et al. (2005) (§4). In addition to deriving 
quasar radio properties as a function of redshift and optical 
color, we reexamine the issue of whether the radio-loudness 
distribution is bimodal. We also explore the distinction be- 
tween broad absorption line (BAL) and non-BAL objects, 
finding the surprising result that BAL quasars have a higher 
mean flux density and radio loudness than non-BAL objects 
below 2 mJy. We conclude (§5) with a summary of the impli- 
cations of our results, and preview other applications of our 
stacking procedure. 

2. MEAN VERSUS MEDIAN STACKING PROCEDURES 

We have explored two different methods for stacking sub- 
threshold FIRST images, one using the mean of each pixel in 
the stack and the other using the median. Both approaches 
have advantages and disadvantages. The mean flux density in 
a stacked image is mathematically simple and is easily inter- 
pretable. However, it is very sensitive to the rare outliers in 
the distribution. The presence of a bright source in the stack, 
either at the image center or in the periphery, makes itself 
obvious in the summed image. Moreover, noise outliers can 
also cause problems, as a minority of very noisy images may 
substantially raise the noise in the mean image. The outlier 



problem can be addressed by testing each image in the stack, 
discarding sources that are actually above the FIRST detec- 
tion threshold and/or discarding images that exceed some rms 
noise threshold. However, the resulting mean is sensitive to 
the exact value of the discard threshholds and hence does not 
provide a very robust measurement. 

The alternative is to determine the median value of the 
stacked images. The obvious advantage is the insensitivity 
of the median to outliers, since the median is well known to 
be robust in the presence of non-Gaussian distributions (e.g., 
Gott et al. 2001). Therefore all of the data can be utilized, 
eliminating the need to impose an arbitrary cutoff to the dis- 
tribution. However, the interpretation of the median value for 
low signal-to-noise (S/N) data is not straightforward. For high 
S/N data, the median is simply the value at the midpoint of 
the distribution. But in the case of low S/N data, the value 
obtained by taking the median is shifted from the true me- 
dian toward the 'local' mean value. The degree of the shift 
depends on the amplitude of the noise; as the noise increases, 
the recovered value approaches the 'local' mean, where the 
'local' mean is the mean of the values within approximately 
one rms of the median. Hence the recovered median value 
can depend on both the intrinsic distribution of the parame- 
ter and the noise level. Some concrete examples may help 
to illuminate this effect. Figure ^ shows the dependence of 
the measured median value (computed using numerical inte- 
gration) on the noise level for two asymmetrical distributions. 
First consider a simple distribution consisting of two Gaus- 
sians centered at xi =0 and X2 = 1 with widths o\ =0.1 and 
(72 = 0.5 (Fig.[^). The Gaussians are normalized to have equal 
integral amplitudes, but because the first is much narrower, its 
peak is higher by a factor of o^/di = 5. The mean of the dis- 
tribution is midway between the Gaussians at (x) = 0.5, but 
the median is dominated by the much better-localized nar- 
row component and falls at median(;t:) = 0.167. If samples 
are drawn from this distribution with additive Gaussian mea- 
surement noise, the peaks of both distributions get broadened. 
In the limit where the noise is much larger than X2 —Xi, the 
median value converges to the mean (x). The transition with 
increasing values of the noise rms is shown in Figure^c). 

The general shape of the transition in Figure^c) is typical 
for simple skewed distributions (e.g., power laws, exponen- 
tials, etc.). However, more complicated distributions display 
a more complex dependence on the noise level. Figure ^b) 
shows a distribution composed of two one-sided exponen- 
tials, P(x)dx = exp(-x / h) / h, with x > 0. The first exponen- 
tial drops rapidly, with a scale height hi = 1, while the second 
drops much more slowly, hq. = 1000. The two components are 
normalized so that the first contains the vast majority of the 
sources, with the integrated amplitude of the second compo- 
nent being only 0.05% of the total. When Gaussian noise is 
added, there are three separate regimes of behavior (Fig.^T). 
When the rms noise is much smaller than either exponen- 
tial scale, the true median is recovered: median(x) = 0.694, 
just slightly above the median computed for component 1 
only (ln2 = 0.693). For very large rms noise levels the mea- 
sured median converges to the mean for the whole distribution 
((x) = 1.50). But for intermediate values of the rms around 
unity, there is an inflection where the measured median value 
pauses at a value of median(x) ~ 1. This is explained by 
the fact that the dominant exponential with hi = 1 is itself a 
skewed distribution having a mean (x) i = 1 . 

It is worth noting that the standard arithmetic mean is 
also of limited utility in the presence of complex multiple- 
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FIG. 1 . — The effect of measurement noise on the measured median value for a skewed distribution, (a) A distribution consisting of two equal Gaussians, a 
wide component at x = 1 and a narrow component at x = 0. (b) A distribution composed of two exponentials having very different scale heights with the broad 
component including only a small fraction of the population. The probability has been multiplied by x for a better visual display of the distribution plotted versus 
log x. (c) The median for the double Gaussian distribution with noise added shows a smooth transition. When the noise is small, we recover the true median, and 
as the rms noise becomes comparable to the separation of the components, the value converges to the mean for the distribution, (d) The median for the double 
exponential distribution with noise presents a transition from the true median to the mean that pauses at the mean for the dominant (narrower) population. 



component, strong-tailed distributions like that in Figure^b). 
Even for those distributions the median is generally a better 
match to one's intuitive concept of the "typical" value of the 
distribution. 

Despite these complications (about which we have found 
little discussion in the astronomical literature), we believe that 
the median is distinctly preferable to the mean for stacking 
our FIRST survey images. In our tests, the robust median 
calculation produces significantly more stable results with 
lower noise, while giving very similar measured values for 
the fluxes. We are in a limit where almost all the values in our 
sample are small compared with the noise, so it is straight- 
forward to interpret our median stack measurements as repre- 
sentative of the mean for the population of sources with flux 
densities fainter than a few times the FIRST rms (i.e., a few 
x 0.145 mJy). Throughout the paper, we refer to the median- 
derived approximate mean interchangeably as the "median" 
or "average" of the quantity of interest. 



3. CALIBRATION OF THE STACKING PROCEDURE 

3.1. Introduction 

As an aperture synthesis interferometer, the Very Large Ar- 
ray 7 samples the Fourier transform of the radio brightness dis- 
tribution on the sky. To obtain a sky image requires trans- 
forming to the image plane with incomplete information. The 
165-second snapshots that comprise the FIRST survey are 
particularly problematic in this regard: we typically obtain 
~ 30,000 visibility points and transform them into images 
with ~5x 10 5 resolution elements. The nonlinear algorithm 
'CLEAN' (Hogbom 1974; Clark 1980) is used to minimize 
artifacts such as the diffraction spikes produced by the VLA 
geometry and the grating rings imposed by the minimum an- 

7 The Very Large Array is an instrument of the National Radio Astronomy 
Observatory, a facility of the National Science Foundation operated under 
cooperative agreement by Associated Universities, Inc. 
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tenna spacings employed. 

One consequence of this process, discovered in the course 
of the NVSS (Condon et al. 1998) and FIRST surveys, is 
"CLEAN bias". This incompletely understood phenomenon 
steals flux from the above-threshold sources and redistributes 
it around the field. The magnitude of the bias is dependent 
on the rms noise in the image (it increases as noise increases), 
the off-axis angle (it decreases in consort with the primary 
beam pattern), and the source extent (extended sources lose 
more flux). The FIRST and NVSS surveys took considerable 
pains to calibrate CLEAN bias, concluding, respectively, that 
it had values of 0.25 mJy/beam and 0.30 mJy/beam. (It is 
unsurprising that the different resolutions, integration times, 
and analysis procedures of the two surveys produced slightly 
different results). 

While the sources of interest in a stacking analysis are 
sub-threshold and therefore, by definition, have not been 
CLEANed, we have taken the discovery of CLEAN bias as 
a cautionary tale and have examined in detail the behavior of 
our images subjected to the stacking process. We find that we 
do not recover the full flux density of either artificial sources 
inserted into the images or real sub-threshold sources. We de- 
scribe here our calibration of this phenomenon, which we dub 
"snapshot bias". 

3.2. Artificial source tests 

The Astronomical Image Processing System (AIPS) used to 
reduce VLA data includes a task UVMOD that allows the user 
to insert artificial sources into a uv database in order to test 
the fidelity of the analysis. Since in this case we are interested 
in sources far below the detection threshold, a large number 
of artificial sources is required. The median rms in coadded 
FIRST fields is 145 /iJy; thus, to achieve an uncertainty of 
~ 10% in the measured flux density of, say, 40 /iJy sources 
requires the addition of more than 1200 individual sources. 

For our initial attempt to measure the bias, we inserted one 
hundred 40 /iJy sources placed in a regular square grid into 
each of 100 FIRST fields. This approach allowed us to mini- 
mize the number of maps we needed to make. However, this 
failed as a consequence of the interference of the sidelobe pat- 
terns that even these very faint sources produce. Constructive 
and destructive interference of sidelobes led to noise in the 
stacked images that varied strongly and systematically at dif- 
ferent artificial source grid positions. We concluded that this 
test was not sufficient to measure the bias for random source 
locations. 

Our ultimate artificial source test involved placing four 
40 jiiy sources in each of 400 FIRST uv datasets. The sources 
were placed at the corners of a square of side ~ 2' centered on 
the image. We then CLEANed the 400 images, extracted 1' 
cutouts around each of the fake source locations, and stacked 
the cutouts to find the median flux density. Since artificial 
source locations were not screened in advance, they ocassion- 
ally fell on or near the location of a real radio source; the 
median algorithm effectively rejected the contaminated pixels 
in those cases (§2). Source parameters were derived by fit- 
ting an elliptical Gaussian to the stacked image as is done for 
source extraction in the real images. To improve the quality 
of the fits, regions around the diffraction spikes in the VLA 
dirty beam (see Fig.|4]below for an example) are masked out. 
The process was repeated for artificial sources with a peak 
flux density of 80 /iJy. 

The results are presented in Table 1 . The recovered median 
peak flux densities for the 40 /iJy and 80 /iJy sources were 



TABLE 1 

Flux Density Bias in Stacked FIRST Images 



True Flux 3 


Stack Median b 


Bias c 


No. Images' 1 


GuJy) 










Artificial Inserted Sources 


40 


35 


5±3 


1600 


80 


61 


19 ±4 


1600 


First-Look Survey 


182 


134 


48 ± 15 


144 


198 


144 


54 ± 14 


145 


219 


154 


65 ± 12 


144 


243 


140 


104 ± 15 


145 


277 


203 


74 ± 10 


144 


320 


231 


89 ±27 


145 


385 


288 


98 ± 11 


144 


492 


313 


180 ± 14 


145 


733 


503 


230 ± 22 


144 


1300 c 


1047 


253 ± 41 


145 


COSMOS Survey 


200 


192 


8 ±20 


37 


328 


234 


93 ±61 


38 


594 


354 


240 ± 12 


37 


1086 c 


893 


193 ± 133 


38 



"Mean peak flux density for sources in flux bin. 
b Median peak flux density for FIRST image stack. 
c Snapshot bias (underestimate of true flux) and rms uncertainty. 
d Number of sources and images in this bin. 

e This value is the median instead of the mean because the noise in the indi- 
vidual images is small compared with the bin's flux density range. 



33 /iJy and 60 /iJy, respectively. The persistence of missing 
flux reminiscent of the CLEAN bias at flux densities far be- 
low those that experience CLEANing is a surprise; as shown 
below, however, this result is confirmed by stacking results on 
faint radio sources derived from deep, full-synthesis images. 

3.3. Recovering real sub-threshold sources 

An alternative to using artificial sources is to stack real ra- 
dio sources detected in very deep VLA surveys. The First- 
Look Survey (FLS - Condon et al. 2003) covered the 5 deg 2 
of the Spitzer First-Look fields using the VLA B configura- 
tion; it achieved a mean rms of 23 /iJy beam" 1 and detected 
3565 sources down to a flux density of 115 fiJy. We ran our 
FIRST survey source extraction routine HAPPY (White et al. 
1997) on the publicly available FLS radio images and con- 
structed a catalog of 1445 point-like sources (deconvolved 
source size < 2.5" with the 5.0" beam) with flux densities 
ranging from 0.17 to 3.0 mJy; we chose a higher (7er) source 
detection threshold to minimize the uncertainties on the indi- 
vidual source flux densities. We grouped the sources into ten 
equally populated flux density bins with mean 8 FLS fluxes 
ranging from 182 /iJy to 1300 /iJy. We then extracted 1' 
cutouts around each of these sources in the FIRST images and 
compared the true mean flux density (from our FLS catalog) 
with the median stacked flux density in each bin. Source pa- 
rameters were derived by fitting an elliptical Gaussian to the 

8 The mean was used for the sub-threshold sources, since this is the value 
to which our median stacking converges (see §2), but the median was used 
for the final bin, which contains detected FIRST sources. 



Stacked FIRST Quasars 



5 




0.05 - 



o.oo r* , , , ,_j 

0.1 1.0 
True Flux Density [mjy] 

FIG. 2. — Snapshot bias for stacked FIRST images as a function of the 
true flux density. The bias is the difference between the true flux density 
and the flux density measured in the median image. The sample is limited to 
point sources (FWHM < 2.5"), and fitted peak values are used for the flux 
densities. Data are shown for both the The First-Look Survey (Condon et al. 
2003; circles) and artificial 40 and 80 /jjy sources (diamonds). The vertical 
bars indicate la errors on the stacked flux density, and the horizontal bars 
represent the range of flux densities for each bin. The solid line is a lineal' 
model in which the bias is 29% of the true flux density. The dashed line is 
the 0.25 mjy CLEAN bias measured for sources bright enough to be in the 
FIRST catalog. 



stacked image as described above. 

We performed a similar analysis using data from the COS- 
MOS survey's pilot program (Schinnerer et al. 2004), com- 
prised of seven VLA pointings in the A configuration that 
reached rms values ranging from 36 to 46 p,Jy. The results 
are consistent with the FLS survey, but the uncertainties are 
much larger because the COSMOS sample has only one tenth 
as many sources as the FLS sample. 

The results are displayed in Table 1 and plotted in Fig- 
ure 13 For sources brighter than 0.75 mjy, the mean deficit 
is consistent with the 0.25 mjy CLEAN bias we have added 
to all above-threshold FIRST sources (see above). Below 
0.75 mjy, however, the flux deficit changes character, and 
is well-represented by a constant fractional offset, with the 
stacked image yielding a value 71% that of the true mean flux 
density for each bin. 

It is perhaps unsurprising that the bias should change near 
the 0.75 mjy threshold that divides brighter sources that were 
CLEANed during the FIRST image processing from fainter 
sources that have not been CLEANed. The continuity be- 
tween the bias for sub-threshold sources and that for super- 
threshold sources suggests that they are aspects of the same 
phenonmenon. A connection between the bias and the depth 
of CLEANing is well established for brighter sources (Becker, 
White & Helfand 1995; Condon et al. 1998), so we think it 
likely (but by no means certain) that the snapshot bias is also 
created by the non-linear CLEAN process. 

However, we are clueless as to why the relationship has the 
particular form we observe. According to Cornwell, Braun & 
Briggs (1999), "to date no one has succeeded in producing a 



noise analysis of CLEAN itself", so we are not alone in being 
mystified. While we cannot offer a theoretical explanation for 
snapshot bias, we will use the simple empirical bias correc- 
tion: 

5 p , corr =min (1.40 S P ,S P +0.25 mJy) , (1) 

where S p is the fitted peak flux density measured from the 
median stack. 

Although Eq. Q was derived from elliptical Gaussian fits 
to the stacked images, we find it applies equally well if the 
brightness of the central pixel in the median image is used to 
estimate the peak flux density. We use both approaches below 
in the analysis of the quasar sample. 

Note that it is most fortunate that the bias is a constant frac- 
tion of the flux, since that means that it can be corrected in 
the stacked image. That would not be true if, for example, it 
were a quadratic function of flux, since in that case the bias 
in the summed image would depend on the detailed distri- 
bution of contributing fluxes (which is unknown). But since 
all faint sources have the same bias correction multiplier, the 
bias correction can be appropriately applied to the stacked im- 
age instead of the individual images. In fact, it can be ap- 
plied pixel-by-pixel to the stacked image by simply multiply- 
ing each pixel in the image by 1 .40. 

3.4. Stacking White Dwarfs, a Radio-Silent Population 

In the remainder of the paper we discuss the results of stack- 
ing quasars divided into groups using many different param- 
eters (redshift, optical luminosity, etc.) In all of the quasar 
subpopulations we stack, we always detect a positive signal. 
In order to allay concern that our algorithm somehow guaran- 
tees a detection, we have stacked 2,412 white dwarfs from the 
SDSS DR1 white dwarf catalog (Kleinman et al. 2004). As 
expected, the stacked image shows no hint of any source. The 
image rms is 3.6 /iJy, comparable to the value expected from 
the typical FIRST rms of 145 ^Jy divided by VN = 49. 

4. THE RADIO PROPERTIES OF UNDETECTED QUASARS 

Although radio emission was the defining feature of the first 
quasars, more than four decades of effort has failed to es- 
tablish predictive models for quasar radio properties. While 
~ 10% of quasars are relatively bright at centimeter wave- 
lengths (> 1 mJy) and thus are readily detected, the radio 
emission from most quasars falls well below the limits of all 
large-area radio surveys. Even deeper surveys that cover sev- 
eral square degrees of sky only detect ~ 50% of quasars. For 
example, by examining the images from the FLS radio sur- 
vey described above (Condon et al. 2003), we detect 36 of 
72 SDSS quasars to a limiting flux density of ~ 0.09 mJy. 
Figure [5] shows the fraction of detected sources as a function 
of flux density. The width of the shaded bands represent the 
lcr and 90% confidence uncertainties; it is apparent that even 
this, the largest of all radio surveys to this depth, is inadequate 
for determining accurately the detected fraction as a function 
of flux density, let alone for understanding how radio emis- 
sion depends on redshift, absolute magnitude, the presence of 
broad absorption lines, etc. For the forseeable future, 90% 
of quasars will remain undetected at radio wavelengths. By 
using image stacking with the FIRST survey, however, one 
can begin to quantify the statistical properties of quasar radio 
emission at all flux density levels. 

4. 1 . The radio properties of SDSS quasars 
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FIG. 3. — The cumulative fraction of radio-detected quasars as a func- 
tion of 20 cm flux density. The shaded bands represent the ±1<t and ±90% 
uncertainties derived using the FLS images for the 72 SDSS quasars in the 
FLS survey area. The cross at 0.75 mJy represents the fraction of all SDSS 
quasars detected in the FIRST survey above this flux density, while the dot 
at 50% fraction indicates the value of the median flux density derived from 
our stacking analysis. The general agreement of the latter with the fraction 
of directly detected quasars at these flux densities offers validation of our ap- 
proach, but the noise (which is comparable to the symbol size) is far smaller 
in the values derived from the stacked images. 



As a starting point we use the largest existing quasar survey 
as reported in the SDSS DR3 catalog (Schneider et al. 2005), 
which contains 46,420 spectroscopically identified quasars. 
Of these 41,295 fall in regions covered by the FIRST survey. 
Constructing a median stack of the entire sample yields the 
image shown as the inset in Figure |4] This high signal-to- 
noise (~ 75 : 1) image shows a compact source centered on 
the nominal quasar(s') position; a two-dimensional Gaussian 
fit yields a peak raw flux density of 80 /iJy, or roughly 50% 
of the rms of an individual FIRST image. Multiplying this by 
1.40 to correct for the snapshot bias (Eq. 1) gives a peak flux 
density of 112 ± 1.5 /iJy. The fluxes plotted in Fig.0]and all 
fluxes quoted hereafter have been corrected for the bias. 

The six positive-flux radial spokes and interspersed nega- 
tive features are characteristic of the VLA sidelobe pattern; 
since the vast majority (~ 93%) of sources contributing to this 
image are below the FIRST detection threshold and are there- 
fore not CLEANed, this sidelobe pattern is expected. Note 
that nearly all FIRST fields are observed near the meridian 
so that sidelobes from different fields align well. The pixel- 
by-pixel radial profile in Figure0]shows a FWHM of ~ 7.0", 
slightly larger than the size expected for a point source ob- 
served in the VLA B configuration at 20 cm. 9 The shaded 
horizontal band indicates the ±ler values derived for median 
statistics (Gott et al. 2001). The intrinsic radio source size 
implied by the extended emission is affected by the VLA PSF 

9 The FIRST survey's cataloged sources have a FWHM of 5 .4" as a con- 
sequence of the fact that the CLEANed FIRST images are convolved with a 
CLEAN beam with that value; this is typically slightly larger that the dirty 
beam size to accomodate images observed away from the zenith where the 
synthesized beam shape is larger than the nominal B-configuration value. 



size, which depends on the distance of the source from the 
zenith. Since most FIRST fields were observed close to the 
meridian, the zenith distance is a simple function of the source 
declination. Figure[3]shows the measured radio sizes in image 
stacks separated into nine zenith-distance bins. The increase 
in size at high zenith distances is explained by the increase in 
the VLA beam size. Since the many quasars being averaged 
for this measurement are randomly oriented, the stacked radio 
image is expected to be symmetrical, and asymmetries are ex- 
plained by beam effects. Both the distribution with zenith dis- 
tance and the fitted size for the image in Figure |4] (6 "4 x 7"0 
FWHM) are consistent with a symmetrical quasar image hav- 
ing a mean source size of 3 "5. This is a bit larger than the size 
of quasars at the 1 mJy detection limit of the FIRST survey. 
Fitting the mean stack for the 679 quasars with central flux 
densities between 1 and 2 mJy yields a size of 5"8 x 6"2, im- 
plying an underlying source size of 2."0 x 3"0 when the beam 
size is deconvolved. 

The median flux density of ~ 110 /iJy is reasonably con- 
sistent with that found for the directly detected quasars within 
the FLS sample (Fig. [3}, though it is slightly higher than the 
value in the FLS field (74 ± 20 ^Jy). It is likely that this differ- 
ence is mainly the result of sample variance in the FLS field. 
If we stack the FIRST images for only the quasars in the FLS 
fields, the flux density is 76 ± 26 mJy, in good agreement with 
the measurement from the FLS images. 

4.2. Variation of radio properties with optical luminosity 

Dividing the SDSS quasars into ten redshift bins, we see 
that the median flux density declines monotonically to z = 2 
(Fig. [6). At z = 2.25 there is a noticeable jump in the radio 
flux, which is a consequence of a confluence of effects driven 
by the sharply declining efficiency of the SDSS quasar se- 
lection algorithm (because the colors of z = 2-3 quasars are 
similar to stars; Richards et al. 2001, 2002) combined with an 
interesting dependence of th e rad io emission on optical color 
(discussed further below in 34.3b . 

The decline of flux with redshift is slower than the expected 
scaling as the inverse of the luminosity-distance squared be- 
cause the SDSS sample is flux-limited and so detects increas- 
ingly luminous objects as the redshift increases. An inter- 
esting point is that although the FIRST catalog is also flux- 
limited, the stacked FIRST data are not flux-limited. All 
sources get included in the stack regardless of their radio 
brightnesses. Consequently these data do not suffer from the 
usual bias against faint sources in the radio; only the optical 
flux limit introduces such a bias. The radio luminosities are 
biased toward brighter values at high redshifts only insofar as 
the radio and optical luminosities are correlated. 

The correlation between radio and optical luminosities does 
introduce complications in interpreting our results. Figure 
displays the radio flux as a function of SDSS /-band magni- 
tude. Optically bright sources are far more likely to be radio 
bright; in fact, for the brightest quasars with i < 16, the me- 
dian radio flux density approaches the 1 mJy detection limit 
for the FIRST survey. This is consistent with the conclusion 
from the FIRST Bright Quasar Survey that FIRST detects 
most V ~ 15 quasars (White et al. 2000.) But the potential 
entanglement of redshift, absolute magnitude, and evolution 
makes it difficult to understand the physical implications of 
this correlation. 

We have concluded that the best approach is to correct for 
the correlation between absolute magnitude and radio lumi- 
nosity before attempting to understand the variation in ra- 
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dio brightness with secondary parameters. We compute the 
2500 A absolute magnitude, Muv> by applying a redshift- 
dependent /^-correction derived using the Vanden Berk et 
al. (2001) composite SDSS quasar spectrum (Fig. [8). The 
/T-correction is applied to the filter closest to 2500 A at 
the quasar redshift. The observed 20 cm flux densities are 
converted to a rest-frame 5 GHz (6 cm) radio luminosity, 
Lr(5GHz), using a spectral index of a = -0.5. The redshift 
is converted to luminosity distance using a standard WMAP 
cosmology (fi m = 0.3, A = 0.7, h = 0.7). The particular choice 
of rest frame wavelengths facilitates comparison of our results 
with previous studies using the radio-loudness parameter/?*, 
defined by Stocke et al. (1992) as the ratio of the 2500 A and 
5 GHz flux densities. We convert each FIRST cutout image 
to radio luminosity units using the known quasar redshift and 
then stack those scaled images to compute median radio lu- 
minosities. Figure|9]shows the very close correlation between 
Muv and Lr, which is well fitted by a power-law: 

logL fi = 0.54-0.339(M £/v +25) , (2) 

where Lr is in units of 10 30 ergss _1 Hz _1 . If the radio lumi- 



nosity were simply proportional to the optical luminosity, the 
slope would be steeper (-0.4 instead of -0.339). This slope 
implies Lr ~ ZM 5 ; the radio loudness R* is a declining func- 
tion of optical luminosity, with the most luminous sources 
(Muv = -28.5) having R* values that are lower by a factor 
of 3 compared with the least luminous sources (Mj/y = -20). 
The absolute magnitude is strongly correlated with redshift, 
but if the sample is divided into redshift intervals we find that 
a similar Lr versus Muv correlation applies at all redshifts 
(Fig.llOL For this test we have restricted the sample of DR3 
quasars to those selected using the primary or high-z target- 
ing criteria (Schneider et al. 2005). The primary selection 
used ugri colors to identify quasar candidates at z < 3 with 
magnitudes i < 19.1; it includes 25,511 objects in regions 
covered by FIRST. The high-z criterion used griz colors to 
identify candidates to fainter levels (z < 20.2) and at redshifts 
greater than 3; our sample includes 2,412 such objects. The 
bulk of the remaining DR3 quasars were selected using var- 
ious serendipity criteria. We exclude them here because ob- 
jects so selected have unusual radio properties (as discussed 
further below in §4.31 ) 
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In order to remove this strong radio-optical correlation, we 
scale the radio properties to the reference absolute magnitude 
Muv = _ 25. This is accomplished simply by multiplying the 
FIRST cutout by an Ma y -dependent factor: 

logSW = 0.339(M c/v +25)+logS , (3) 

where S is the original radio flux density. The adjustment to 
the radio luminosity, Lr.m, is similar, and the adjusted radio- 
loudness ratio is: 

log^ = -0.061(Mt/ V + 25) + log/?* . (4) 

In all cases the M subscript indicates that the quantity has been 
adjusted for the absolute magnitude dependence. 

4.3. Variation of radio properties with redshift and color 
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FIG. 8. — iT-correction as a function of redshift to convert observed 
SDSS magnitudes to the magnitude at 2500 A rest wavelength, derived us- 
ing the SDSS composite quasar spectrum (Vanden Berk et al. 2001). The 
^-correction is added to the magnitude in the SDSS filter closest to 2500 A 
(shown by the vertical bands). The redshift-dependent bandwidth-stretching 
factor, — 2.51og(l +z), has been omitted here to make the plot clearer. 



Figure II ll a) shows the redshift dependence of the radio 
loudness R* M after adjusing for the dependence on absolute 
magnitude. There is only modest evolution in this quantity, 
with the radio-loudness being a factor of two higher at z = 5 
than at z = 0. The radio properties of typical quasars have 
changed little since the universe was one billion years old. 

The picture changes, however, if we separate the SDSS 
DR3 sample according to the criteria used to select the can- 
didate quasars for spectroscopic observations (Fig. II lb). The 
quasars selected using the high-z criterion, which are redder 
and fainter than the primary candidates, are brighter in the ra- 
dio. This is most noticeable at low redshifts (z < 2), where the 
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difference in brightness is a factor of 4. Even at high redshifts 
(z > 3) a slight difference persists; at least part of the slow rise 
in R'm toward higher redshifts (Fig. II Ih) appears to be created 
by the transition in the SDSS sample frompn'raary-dominated 
selection for z < 3 to /z/g/j-z-dominated selection for z > 3. 

Quasars selected using other criteria (serendipity, ROSAT, 
FIRST, stars, etc., as described in Schneider et al. 2005) are 
also systematically radio louder. One might be tempted to 
ascribe this to the use of the FIRST catalog in selecting some 
of these candidates; however, that introduces at most a very 
minor bias toward higher R^ values. Only 279 of the 13,372 
sources selected using other criteria are flagged in the DR3 
catalog as FIRST sources, and excluding them reduces the 
median radio loudness only slightly from log/?^ = 0.11 to 
0.08. (The robustness of the median to the presence of a rare 
admixture of bright sources is of course the reason we choose 
to use it.) Similarly, excluding ROSAT-selected sources — 
since radio emission is known to be more common among X- 
ray quasars (e.g., Green et al. 1995) — also leads to only a 
very small reduction in R^- We conclude that there must be 
another explanation for the different radio properties of the 
variously selected quasar samples. 

One possible contributor is the anti-correlation between ra- 
dio loudness and apparent magnitude (Fig. 1121 . The optically 
fainter sources are radio-louder, even after the Muv adjust- 
ment. Quasar candidates selected using the primary criterion 
are on average 1 magnitude brighter than those selected using 
other criteria (i = 18.6 versus 19.6). But that accounts for a 
difference in log R^ of only 0. 15 and so does not explain the 
bulk of the difference between the samples. 

The most important underlying correlation that leads to dif- 
ferences between the different SDSS samples is a strong de- 
pendence of radio loudness on optical color (Fig. ll3l >. Since 
quasar colors change systematically with redshift as various 
emission lines move through the SDSS filters, we have sub- 
tracted the color of the SDSS composite quasar spectrum 
(Vanden Berk et al. 2001) from the observed colors. That re- 



duces the scatter in colors and makes the expected color zero 
for a quasar that resembles the composite. We see a strik- 
ing correlation: quasars that are either bluer or redder than 
the standard color are brighter in the radio, and substantially 
redder objects (with g — r > 0.8 magnitudes) are brighter by a 
factor of ~ 10 than quasars with typical colors. The tendency 
of radio-loud quasars to have a larger scatter in their optical 
colors has been noted before (Richards et al. 2001, Ivezic et 
al. 2002), although it has never been so clearly seen as in this 
analysis. Not only are the red quasars radio-louder, but their 
median flux densities are also far higher (Fig. I14l >. The in- 
crease in Rlf for red objects is due primarily to brighter radio 
fluxes, not to fainter optical magnitudes (which might also be 
expected if the reddening is due to dust extinction.) Of the 
factor of 10 variation seen in R^, a factor of 4 is attributable 
to brighter radio flux densities and a factor of 2.5 to fainter op- 
tical fluxes. Note that the reddest sources have median radio 
flux densities of nearly 0.4 mJy, tantalizingly close to detec- 
tion by the FIRST survey. 

Figure^] also shows the distribution of colors for quasars 
selected using the various SDSS candidate criteria. Quasars 
selected using the primary criterion are much more concen- 
trated toward the typical (zero) colors than are objects selected 
by either the high-z or other criteria. This makes a signif- 
icant contribution to the radio-loudness differences between 
the various samples (Fig. II lb). The color differences between 
the primary and high-z samples, when folded through the cor- 
relation in Figure[Dl lead to a difference in log R* M of ~ 0.6 
between the samples for low-redshift quasars (z < 1.5). That 
accounts for most of the difference between the samples seen 
in Figure fTTT b). 

The radio-color correlation also accounts for the bump seen 
in the flux density at z= 2.25 (Fig . |6j . The efficiency of the 
SDSS color selection for quasars declines sharply in the red- 
shift range 2.4 < z < 3 because the locus of normal quasar col- 
ors crosses the stellar locus in the SDSS color space (Richards 
et al. 2001, 2002). The SDSS DR3 catalog contains far fewer 
objects in this redshift range than might be expected based on 
the sensitivity of the survey. Moreover, the quasars that are 
included in the catalog are dominated by objects with unusual 
colors compared with the composite spectrum, since such ob- 
jects do not resemble stars and so can be selected by the usual 
SDSS criteria. The jump in the radio flux over this redshift 
range is created by the selection of a larger fraction of redder 
quasars that have brighter radio emission than normal quasars. 

The observed increase in radio emission for bluer quasars 
can be understood in the context of the unified model for AGN 
(Urry & Padovani 1995) as being due to a beamed blazar com- 
ponent affecting both the optical and radio emission. The 
brightening for redder sources could also be attributed to red 
optical synchrotron emission, but might instead be explained 
as an evolutionary effect where dusty quasars are more likely 
to be low-level radio emitters. This is likely to be a very useful 
clue to further understanding of the origins of radio emission 
in AGN. 

4.4. The radio-loudness dichotomy 

Our absolute magnitude-adjusted radio loudness parameter 
R'm can be used to re-examine the issue of whether the radio 
loudness distribution is bimodal. While all observers agree 
that there is a highly non-Gaussian tail toward high R* val- 
ues, White et al. (2000) and Cirasuolo et al. (2003a,b) did not 
see evidence for a truly bimodal distribution with two peaks. 
Ivezic et al. (2002) did find a secondary peak, though their 
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FIG. 10. — The radio luminosity as a function of absolute magnitude with the quasar sample divided into redshift intervals. Only DR3 quasars selected using 
the primary or high-z selection criteria are included. The upper left panel shows the combined sample; the number in parentheses gives the number of sources in 
each redshift interval. Despite the strong absolute magnitude-redshift correlation, which can make it difficult to separate dependencies on the two variables, it is 
clear that the quasars in all the redshift bins follow the same Lr versus Muv relationship with at most modest variations. 
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methodology was questioned by Cirasuolo et al. Ivezic et 
al. (2004) subsequently applied the Cirasuolo et al. (2003b) 
approach to a large sample of SDSS quasar candidates and 
claimed conclusive evidence for a double-peaked distribution. 

Figure^]shows our distribution for the radio-loudness pa- 
rameter. Note that the R'^ parameter includes both redshift- 
dependent ^-corrections (as recommended by Ivezic et al. 
2004) and our absolute-magnitude adjustment. The overall 
distribution (Fig. I15h ) clearly does show a secondary peak, 
although the contrast in the valley between the peaks is con- 
siderably less than the factor of two found by Ivezic et al. 
(2004). The peak is also at a considerably lower log R* value 
after the absolute magnitude adjustment (logR^ ~ 1.55 in- 
stead of log R* ~ 1 .9). 

An exploration of the dependence of the R^ distribution 
on other parameters reveals a complex situation. The radio- 
loud tail is considerable weaker at low redshifts (z < 0.5; 
Fig. 115b) but is especially strong in the intermediate redshift 
range (2.5 < z < 3; Fig.ll5t) where the color-section effects 
discussed above are dominant. For red sources, the valley dis- 
appears altogether with the resulting distribution being shifted 
by a factor of ~ 3 toward higher radio-loudness (Fig. I15t). 
Our conclusion is that there is indeed a double-peaked radio- 
loudness distribution for SDSS DR3 quasars, but that the dis- 
tribution varies dramatically with redshift and color (and other 
parameters). The exact form of the overall distribution is 
likely to have been sculpted by selection effects, which must 
be modeled in detail before the relatively modest 20% dip 
between the radio-loud and radio-quiet quasars can be inter- 
preted. 

4.5. The radio emission ofBAL quasars 

Historically, the strongest claim associating radio properties 
with other quasar attributes was the absence of broad absorp- 
tion lines in radio-loud quasars (Stocke et al. 1992). Becker 
et al. (2000) showed that some radio-loud quasars have BALs, 



although they noted that BALs still appeared to be absent from 
the most extreme radio-loud objects. Recently, Trump et al. 
(2006) released a catalog containing 4784 BAL and near-BAL 
quasars from the SDSS DR3 release. As with quasars in gen- 
eral, most of these objects fall below the detection threshold 
of FIRST, making them an ideal population for stacking stud- 
ies. The FIRST survey covers 4292 of the cataloged BAL 
quasars. 

Traditionally, BAL quasars have been divided into two pri- 
mary subgroups, high ionization (HiBALs) and low ionization 
(LoBALs). The latter are much rarer; in the SDSS/FIRST 
sample there are 3647 HiBALs and only 645 LoBALs. 10 The 
HiBALs are identified primarily on the basis of C IV absorp- 
tion at 1550 A, while the LoBALs are identified primarily 
from Mg II absorption at 2800 A. As a result, the known 
LoBALs tend to be at lower redshifts than the known HiBALs. 
HiBALs can only be identified at redshifts z > 1.7, while 
LoBALs can be recognized at redshifts as low as 0.5. 

In Figure [ToT a), we show the median radio flux density of 
HiBALs, LoBALs, and non-BALs as a function of redshift. 
Interestingly, both classes of BAL quasars are brighter in the 
radio than non-BALs. The absolute-magnitude-adjusted radio 
loudness R^ shows a similar effect (Fig. 116b). indicating that 
this is not due to a difference in the distribution of Mj/y for 
the various classes. The LoBALs are radio-louder by a factor 
2.23 ± 0. 10 (averaged over 0.5 < z < 4) and the HiBALs by a 
factor 1.18 ±0.07 (1.7 < z < 4.3). 

That said, the comparison to the FIRST survey for this 
new large sample of BAL quasars confirms the absence of 
extremely radio-loud BAL QSOs. In Figure WH a), we show 
the cumulative distribution of BAL and non-BAL quasars as 
a function of radio flux density. This plot includes only non- 
BALs with 1 .7 < z < 4, since outside that redshift range the 
absorption lines required for confident identification of non- 
BALs do not fall in the SDSS spectrum window. It is clear 
from the graph that while BALs are not found among the 
brightest radio-emitting quasars, below 2 mJy they are sys- 
tematically brighter than non-BAL objects. 

The disparity remains if we examine the radio-loudness pa- 
rameter instead of the flux density (Fig. 117b). although the 
intermediate brightness HiBALs and non-BALs have R* M dis- 
tributions that are much more similar than their flux distri- 
butions. A two-sided Kolmogorov-Smirnov test shows the 
R'm distributions for HiBALs and non-BALs in Figure fTTl b) 
are different at the 4 x 10~ 5 level of significance. Separate 
tests for the distribution with R'j^ > 14 and the portion with 
R'm < 5.25 show that the bright and faint distributions are both 
discrepant (2 x 10~ 7 and 4 x 10~ 6 , respectively). The differ- 
ences compared with the LoBAL distribution also highly sig- 
nificant (< 4 x 10" 3 ) despite the fact that the LoBAL sample 
is much smaller. The 7?^ distributions for the radio-bright 
(Rlj > 100) HiBAL and LoBAL quasars are statistically in- 
distinguishible. We have also examined the dependence of R^ 
on the BAL quasar catalog's absorption index, which quanti- 
fies the strength and extent of the broad absorption (Trump et 
al. 2006). There is a hint of a decline in R^ at the lowest ab- 
sorption index values, though the size of the effect is modest 
at best (Fig.ll8>. We note in passing that the sample of BAL 
quasars in the SDSS DR3 catalog is strongly biased around 
z ~ 2.5 by selection effects that favor the discovery of objects 

10 The Trump et al. (2006) catalog identifies many sources as both HiBALs 
and LoBALs; we chose to make these categories disjoint by labelling quasars 
as HiBALs only if they are not LoBALs. 
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objects, but the dichotomy is not so clear in different parts of parameter space, (b) Distribution for low-redshift quasars, normalized to show the fraction in each 
bin. The gray line shows the distribution from panel (a) for comparison, (c) Distribution for quasars in the redshift range 2.5 < z < 3, where the SDSS selection 
effects are most severe, (d) Distribution for red quasars. 



with unusual colors. Figure [T^l shows the fraction of BAL 
quasars as a function of redshift. For 2.7 < z < 2.9, almost 
half the objects in the DR3 catalog are BALs! The effect here 
is similar to the bias in favor of the discovery of red quasars 
in this same redshift interval (discussed in jj4.3> . If the broad 
absorption lines change the quasar magnitude in even one of 
the five SDSS filter bands, it is more likely to be recognized as 
having colors inconsistent with the stellar locus. This effect 
has been previously discussed by Hewett & Foltz (2003) and 



Reichard et al. (2003). 

We find the radio dominance of BAL over non-BAL quasars 
difficult to reconcile with claims that BALs are largely the re- 
sult of a preferred orientation (e.g., Murray et al. 1995, Elvis 
2000). In fact, most of the arguments against the orientation 
model to date have been based on radio observations. Zhou 
et al. (2006) argued that radio variability observed in six BAL 
QSOs was strong evidence that at least for some BALs, we 
were looking along the axis of the radio jet. Becker et al. 
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FIG. 16. — Median radio flux density (left) and absolute magnitude-adjusted radio-loudness (right) for HiBAL, LoBAL and non-BAL quasars as a function of 
redshift. Surprisingly, BAL quasars are brighter radio sources than non-BALs, with the effect especially noticeable for low-ionization BALs. 
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FIG. 17. — The cumulative distribution of HiBAL, LoBAL, and non-BAL quasars as a function of flux density limit (left) and magnitude-adjusted radio 
loudness R^ (right). The shaded bands show la uncertainties. The bands converge at fraction 0.5 to the median values derived from our stacking analysis; the 
inset shows an expanded view of that region. There is a deficit of BAL quasars at bright fluxes, but there is an excess at fluxes fainter than ~ 1 .5 mJy for HiBALs 
and ~ 5 mJy for LoBALs. The radio-loudness distribution is similar for non-BALs and HiBALs in the radio-intermediate region (2 < R^ < 10), though a small 
but significant difference remains for radio-quiet (R^ < 2) sources. The LoBALs are much more likely to be radio-intermediate sources than either the HiBALs 
or non-BALs. 



(1997) made similar arguments on the basis of the flat ra- 
dio spectra observed for some BAL quasars. And Gregg et 
al. (2006) used the FR2 radio morphology exhibited by some 
BAL QSOs to argue against the need of a special orientation. 

If the presence of BALs is determined by orientation, the 
greater radio-loudness of BAL QSOs implies that we are look- 



ing closer to, not farther from, the jet axis in quasars with 
BALs. We know of no model that results in higher measured 
radio flux density from the quasar core with greater angular 
distance from the jet direction. Rather, relativistic beaming 
should enhance radio emission at small angles to the quasar 
symmetry axis. This is at odds with conventional orientation 
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FIG. 18. — Median radio-loudness for HiBAL and LoBAL quasars as a function of the BAL absorption index. The top panel shows the distribution of the 
absorption index values for the two types of BALs. The dashed line shows the radio loudness for the non-BAL quasars (which by definition have absorption 
indices of zero.) There is some indication of a drop in radio loudness for low absorption indices, although the effect is not strong. 




FIG . 19. — Fraction of DR3 quasars that are BALs as function of redshift. HiBALs can only be recognized when z > 1.7, which accounts for the large jump at 
that redshift; lower redshift quasars are all of the rare LoBAL type. The spike at z = 2.7 is created by the inefficiency of the SDSS quasar color selection where 
the quasar locus passes close to the stellar locus. 



models that require viewing angles closer to edge-on for BAL 
quasars (e.g., Elvis 2000). 

An alternative explanation is that quasars with low-level ra- 
dio emission in the nucleus have more BAL clouds and con- 
sequently are more likely to show absorption. BAL quasars 
may be in a special evolutionary phase in which low-level ra- 



dio emission confined near the nucleus is accompanied by an 
excess of absorption clouds; when the radio source breaks out 
to become truly radio-loud, it quickly eliminates the clouds 
that are the source of broad absorption lines. This is an evo- 
lutionary unification model (e.g., Lfpari & Terlevich 2006). 
Richards et al. (2004) and Richards (2006) suggest that 
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the sequence from LoBALs to HiBALs to non-BALs may 
be the result of orientation coupled with a gradual transition 
from a radiation-dominated, accretion-disk wind to a MHD- 
dominated wind (Everett et al. 2001; Everett 2005). In this 
view, powerful radio-bright quasars have a tightly collimated 
polar MHD wind in which BAL clouds have a small cover- 
ing factor, while radio-weak objects have more widely dis- 
tributed BAL clouds in an equatorial, radiatively accelerated 
wind coming off the surface of the disk (see Fig. 5 of Richards 
et al. 2004). This is a static unification model in which the 
magnetic field strength (with an associated radio source) acts 
as a second parameter (in addition to orientation). 

This picture does not necessarily predict enhanced low- 
level radio emission from BAL quasars since it nominally pro- 
duces monotonically decreasing radio emission with increas- 
ing BAL strength. However, it could plausibly be modified to 
explain our observations. A weak radio source might elevate 
the disk wind above the equatorial plane, increasing the BAL 
covering factor for lines of sight that are not heavily absorbed 
by the disk, while a stronger radio source would collimate 
the flow and reduce the BAL covering fraction. Any model 
that matches our observations will share this non-monotonic 
behavior in order to produce BALs that are enhanced in the 
presence of low-level radio emission but suppressed by bright 
radio emission. 

5. SUMMARY AND CONCLUSIONS 

We have demonstrated that the average radio properties of 
sources in the FIRST survey area can be derived even for pop- 
ulations in which the individual members have flux densities 
an order of magnitude or more below the typical field rms. 
Our median stacking algorithm is robust, and, following the 
calibration of snapshot bias derived herein, can be used to pro- 
vide quantitative information of average source flux densities 
into the nanoJansky regime. In our application of this algo- 
rithm to the SDSS DR3 quasar catalog, we establish the ra- 
dio properties of quasars as a function of optical luminosity, 
color and redshift. The average radio luminosity correlates 
very well with the optical luminosity, with Lr ~ There 
is a very strong correlation between radio loudness and color, 
with quasars having either bluer or redder colors than the 
norm being brighter in the radio; objects 0.8 magnitudes red- 
der than average in g — r have radio loudness values 10 times 
higher than quasars with typical colors. At faint flux densities, 



BAL quasars actually have higher average radio luminosities 
and radio-loudness parameters than non-BAL objects, a result 
inconsistent with the conventional orientation hypothesis for 
BAL quasars. 

The correlation between radio emission and color is an in- 
triguing clue to the nature of our FIRST-selected red quasars 
(Gregg et al. 2002, White et al. 2003, Glikman et al. 2004a). 
It suggests that a wide-area radio survey only a factor of two 
deeper than the FIRST survey might be capable of detect- 
ing the bulk of the reddened population, which would shed 
light on the still controversial question of what fraction of all 
quasars are highly reddened. 

The success of stacking FIRST images to find the mean ra- 
dio properties of sub-threshold radio sources depends on the 
availability of large target lists. As shown in this paper, the 
SDSS quasar sample is ideal for these purposes. In fact the 
SDSS provides much more than quasars. We are currently 
working on a study of the mean radio properties of SDSS 
narrow-line AGN (deVries et al., in preparation). In that 
paper, we explore the dependence of radio emission on the 
strength of various emission lines, on associated star forma- 
tion, and on black hole mass. We are also examining the radio 
properties of star-forming galaxies taken from the SDSS spec- 
troscopic survey of galaxies (Becker et al., in preparation). 
There is no reason that these studies must be limited to extra- 
galactic samples. Other astrophysical sources of weak r adio 
emission include several classes of stars. (We show in 33.41 
that the average radio flux density of white dwarfs is fainter 
than ~ 10 /iJy.) The number of stars with spectral classifica- 
tion is large enough that it will be feasible to study the radio 
properties as a function of spectral type. 
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